library(dynlm)
## Warning: package 'dynlm' was built under R version 4.0.5
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(zoo)
library(forecast)
## Warning: package 'forecast' was built under R version 4.0.4
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(strucchange)
## Warning: package 'strucchange' was built under R version 4.0.5
## Loading required package: sandwich
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.4     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x stringr::boundary() masks strucchange::boundary()
## x dplyr::filter()     masks stats::filter()
## x dplyr::lag()        masks stats::lag()
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.0.5
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
library(xts)
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(DataCombine)
## Warning: package 'DataCombine' was built under R version 4.0.5
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.0.5
cat("\014")
set.seed(123)

#clear global environment
rm(list = ls())

#below function reads the data from csv file
#Input : String of file
#Output: dataframe
read_data <- function(data){
  data <- read.csv(data)
  #make first column as date
  data[,1] <-as.Date(data[,1],format ="%m-%d-%y")
  data$M2 <- as.numeric(data$M2)
  return(data)
}

df<-read_data("Data.csv")

#convert cpi to inflation and make new dataframe storing in df
df <- change(data = df,Var = "CPI",NewVar = "Inflation",type = "percent")
## 
## Remember to put data in time order before running.
## 
## Lagging CPI by 1 time units.
#convert all columns to monthly time series
inflation <- ts(df$Inflation,frequency=12, start=c(1991,1))
ipilsm <- ts(df$IPILSM,frequency=12, start=c(1991,1))
tbr <- ts(df$TBR,frequency=12, start=c(1991,1))
m2 <- ts(df$M2,frequency=12, start=c(1991,1))
exports <- ts(df$Exports,frequency=12, start=c(1991,1))
imports <- ts(df$Imports,frequency=12, start=c(1991,1))

#split data for train and test and set the last 12 months as a testing partition 
split_inflation <- ts_split(ts.obj = inflation, sample.out = 12)
split_ipilsm <- ts_split(ts.obj = ipilsm, sample.out = 12)
split_tbr <- ts_split(ts.obj = tbr, sample.out = 12)
split_m2 <- ts_split(ts.obj = m2, sample.out = 12)
split_exports <- ts_split(ts.obj = exports, sample.out = 12)
split_imports <- ts_split(ts.obj = imports, sample.out = 12)


#fit auto.arima on each series training data
inflation_fit <- auto.arima(split_inflation[["train"]])
ipilsm_fit <- auto.arima(split_ipilsm[["train"]])
tbr_fit <- auto.arima(split_tbr[["train"]])
m2_fit <- auto.arima(split_m2[["train"]])
exports_fit <- auto.arima(split_exports[["train"]])
imports_fit <- auto.arima(split_imports[["train"]])


#generate predictions for next 12 months
inflation_predictions <- forecast(inflation_fit,h=12,fan = TRUE)
ipilsm_predictions <- forecast(ipilsm_fit,h= 12,fan = TRUE)
tbr_predictions <- forecast(tbr_fit,h=12,fan = TRUE)
m2_predictions <- forecast(m2_fit ,h=12,fan = TRUE)
exports_predictions <- forecast(exports_fit,h=12,fan = TRUE)
imports_predictions <- forecast(imports_fit,h=12,fan = TRUE)

#plot all the forecast with actual data with their fanplots
plot_forecast(inflation_predictions)
plot_forecast(ipilsm_predictions)
plot_forecast(tbr_predictions)
plot_forecast(m2_predictions)
plot_forecast(exports_predictions)
plot_forecast(imports_predictions)
#check residuals of all fitted models and conclude by Ljung-Box test
checkresiduals(inflation_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0)(1,0,2)[12] with non-zero mean
## Q* = 27.369, df = 18, p-value = 0.07232
## 
## Model df: 6.   Total lags used: 24
checkresiduals(ipilsm_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,2)[12] with drift
## Q* = 36.181, df = 18, p-value = 0.006689
## 
## Model df: 6.   Total lags used: 24
checkresiduals(tbr_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)
## Q* = 23.149, df = 24, p-value = 0.511
## 
## Model df: 0.   Total lags used: 24
checkresiduals(m2_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)(0,1,1)[12]
## Q* = 41.622, df = 19, p-value = 0.001993
## 
## Model df: 5.   Total lags used: 24
checkresiduals(exports_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)(0,0,2)[12]
## Q* = 35.443, df = 20, p-value = 0.01787
## 
## Model df: 4.   Total lags used: 24
checkresiduals(imports_fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,0,1)[12]
## Q* = 42.508, df = 22, p-value = 0.005421
## 
## Model df: 2.   Total lags used: 24
#plot training with test data
inflation_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_inflation[["test"]])
## Warning: Removed 1 row(s) containing missing values (geom_path).

ipilsm_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_ipilsm[["test"]])

tbr_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_tbr[["test"]])

m2_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_m2[["test"]])

exports_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_exports[["test"]])

imports_fit %>%
  forecast(h=12) %>%
  autoplot() + autolayer(split_imports[["test"]])